THE LARGE CORE LIMIT OF SPIRAL WAVES IN EXCITABLE 
MEDIA: A NUMERICAL APPROACH 
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Abstract. Wc modify the freezing method introduced by Beyn & Thiimmler, 2004, f° r analyzing 
rigidly rotating spiral waves in excitable media. The proposed method is designed to stably determine 
the rotation frequency and the core radius of rotating spirals, as well as the approximate shape of 
spiral waves in unbounded domains. In particular, we introduce spiral wave boundary conditions 
based on geometric approximations of spiral wave solutions by Archimedean spirals and by involutes 
of circles. We further propose a simple implementation of boundary conditions for the case when the 
inhibitor is non- diffusive, a case which had previously caused spurious oscillations. 

We then utilize the method to numerically analyze the large core limit. The proposed method 
allows us to investigate the case close to criticality where spiral waves acquire infinite core radius r c 
and zero rotation frequency u>, before they begin to develop into retracting fingers. We confirm the 
linear scaling regime of a drift bifurcation for the rotation frequency and the core radius of spiral 
wave solutions close to criticality. This regime is unattainable with conventional numerical methods. 
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1. Introduction. Excitable media are abundant in nature, and appear in phys- 
ical, chemical and biological systems. Prominent examples include cAMP waves in 
slime mold aggregation [39] and intracellular calcium waves [3], electrical waves in 
cardiac and nerve tissue [351 [5], and the auto-catalytic Belousov-Zhabotinsky reac- 
tion |45j . Excitable media involve the interwoven dynamics of so called activators and 
inhibitors. An important class of excitable media contains non-diffusive inhibitors. 
This is the classic situation in cardiac dynamics where the inhibitor consists of (rela- 
tively) immobile ion channels. Excitable media support localized pulses and periodic 
wave trains. In two dimensions rotating vortices (or spirals) are possible |45) . and in 
three dimensions scroll waves occur (35J S3 HH 1301 131] • 

In this work we focus on rigidly rotating spiral waves of single-diffusive excitable me- 
dia. These spiral waves perform a time-periodic motion where the spiral tip moves on 
a perfect circle around an unexcited region, the spiral core. Such spirals are character- 
ized by their rotation frequency u> and their core radius r c . Varying the excitability of 
the medium leads to qualitative changes in the type of motion the spiral is performing. 
For certain parameters rigidly rotating spirals bifurcate via a Hopf bifurcation into 
meandering spirals [3E1 131 122] • A different type of bifurcation occurs at low excitabil- 
ities where rigidly rotating spirals develop into travelling fingers with an overall zero 
curvature. Approaching this bifurcation, the rotation frequency of a rigidly rotating 
spiral wave u decreases down to zero and its core radius r c diverges. Past this bifur- 
cation rigidly rotating spirals do not exist. Whereas the bifurcation to meandering 
spirals is well- understood numerically and theoretically [31 EU HH H3 [33 [3H] , the lat- 
ter one is not. Theoretically, several attempts have been made to study this large core 
limit, using kinematic theory [53J |43j [121 32, 22 , 23l[52l[T3], dynamical systems theory 
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PQ and non-perturbative asymptotics |20j . Computational results, on the other hand, 
are rare. In [20) the large core limit was investigated in terms of the growing velocity 
of the spiral wave tip close to criticality. Using a non-perturbative approach it was 
shown that the growing velocity behaves linearly as a function of the excitability close 
to criticality. This was verified numerically. 

Here we set out to tackle the numerically much more difficult problem of determining 
the rotation frequency u> and the radius r c of the core of a spiral wave in this limit. 
This presents a huge computational challenge. Given current computer power, de- 
termining the rotation frequency and the core radius from direct simulations of the 
underlying reaction-diffusion system is restricted to small core spirals. In the large 
core limit, where the rotation frequency becomes zero and the core radius becomes 
infinite, no reliable results have been obtained so far. 

Our aim here is twofold. First we present a numerical method to study rigidly rotating 
spirals in excitable media which is able to capture the spiral wave even in the large 
core limit without having to solve for computationally expensive large domains. We 
then apply this method to obtain the rotation frequency and the core radius of spiral 
wave solutions in the large core limit, and establish their scaling behaviour close to 
criticality to study the type of bifurcation. 

The numerical method we present is based on the freezing method [H] in which the 
dynamics of an equivariant partial differential equation is split orthogonally into the 
shape dynamics and the dynamics of the associated symmetry group. This method 
has been successfully applied to many types of equivariant systems 1301 SI] • How- 
ever, when applying the freezing method to single-diffusive excitable media several 
problems were encountered [5J HJ H] ■ Firstly, the inhibitor exhibits spurious oscil- 
lations, the amplitude of which increases when approaching the large core limit. It 
was suggested that these oscillations are caused by numerical instabilities linked to 
the mixed hyperbolic-parabolic nature of the problem. The oscillations could be par- 
tially controlled by an upwind-downwind scheme for small core spirals. Secondly, the 
application of Neumann boundary conditions causes the shape of the spiral wave to 
deviate near the boundary from the form expected in an unbounded domain. How- 
ever, it turns out that although the shape is not properly resolved near the boundary, 
the group parameters, i.e. rotation frequency and core radius, are accurately repro- 
duced for small core spirals. This can be understood on a phcnomcnological basis by 
considering that the spiral wave coils shield the solution near the core, and informa- 
tion only flows outwards for rotating spiral wave solutions. 

We will present an implementation of the freezing method which will overcome 
these difficulties. At first, we will identify two separate types of spatial oscillations, 
which are caused by different mechanisms: oscillations in the interior of the compu- 
tational domain linked to numerical instabilities of the non-diffusive inhibitor, and 
oscillations close to the boundary caused by the application of Neumann boundary 
conditions. We eliminate the former by employing a semi-implicit Crank-Nicolson 
scheme and the latter by imposing a subtle implementation of a different type of 
boundary condition. We implement so-called spiral wave boundary conditions, ap- 
proximating the spiral wave by an Archimedean spiral or by an involute of a circle, 
which approximately respect the symmetry of the solution, and produce the correct 
shape of a spiral wave in an effective unbounded domain at the boundary. We find that 
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the spatial oscillations can be controlled if (a) boundary conditions which respect the 
symmetry of the solution arc employed, and (b) the implementation of these bound- 
ary conditions is such that the boundary is coupled to the interior. 

The paper is organized as follows. In Section [2] we present the excitable media model 
under consideration. In Section [3] we describe the original freezing method. This 
method is then modified to suit excitable media with a non-diffusive inhibitor in 
Section 2J where we introduce spiral wave boundary conditions and their discrete 
implementations. In Section [5] we apply our method to study the large core limit 
and determine the scaling behaviour of the rotation frequency and the core radius of 
spiral wave solutions close to criticality for varying excitability e. We conclude with 
a discussion in Section |B1 

2. Model. We consider here the Barkley model [2] for an activator u and an 
inhibitor v described by 

(2.1) d t u = Au + Fiu, v), J-{u, v) — —u(l — u){u — - ) , 

e a 

(2.2) d t v = D v Av + (u - v) . 

Although the numerical method we will describe in Sections [3] and |4] is independent 
of the particular model used, we illustrate some basic properties of excitable media 
with the Barkley model (|2.1j) . Our choice of model is motivated by the fact that 
it incorporates the ingredients of an excitable system in a compact and lucid way. 
Thus, for u s = b/a > the rest state uq = vq = is linearly stable with decay rates 
<Ti = u s /e along the activator direction and a 2 = 1 along the inhibitor direction. 
Perturbing u above the threshold u s (in OD) will lead to growth of u. In the absence 
of the inhibitor v the activator will saturate at u = 1 leading to a bistable system. The 
positive inhibitor growth factor forces the activator to decay back to u = 0. Finally 
also the inhibitor with the refractory time constant 1 will decay back to v = 0. For 
a > b + 1 with b > the system is in OD no longer excitable but instead bistable. 

We have included a diffusion term for the inhibitor v in (|2.1|) . However, we 
will be concentrating on the case of vanishing diffusivity for v with D v = 0. This 
case is more relevant for cardiac dynamics where the inhibitor models the (relatively) 
immobile potassium and sodium ion channels. Moreover, we will see in Section [5] that 
problems of applying the freezing method as proposed in [6] to excitable media are 
caused by the lack of coupling of the boundary and the interior when D v = 0. 

We shall fix in our numerical simulations D v = 0, a = 0.75 and b = 0.01 and vary 
the excitability parameter e if not stated otherwise. 

The Barkley model supports, in a well-defined parameter region [2], rigidly rotating 
spirals. These spiral wave solutions are characterized by their rotation frequency u> 
and their core radius r c . In the large core limit the rotation frequency approaches zero 
and the core radius becomes infinite at a critical value e c of the excitability parameter 
e. At criticality a rigidly rotating spiral wave becomes a finger propagating in the 
transverse direction only, with the speed of the corresponding travelling wave. Finger 
like initial conditions will curl up and eventually develop into rigidly rotating spirals 
below criticality, or into retracting fingers above criticality (see Figure l2~Tj) . Whereas 
the rotation frequency is well defined for rigidly rotating spirals, the definition of the 
core radius is less clear (see [H] for a recent discussion). We define the core of a 
rigidly rotating spiral to be the maximal region which has never been excited during 
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(a) Finger developing into a spiral (b) Retracting finger 

Fig. 2.1. Temporal evolution of a retracting and spiraling finger solution u of the Barkley model 
12. 1\) at 3 different time steps. The motion is clockwise for the spiraling finger and from right to 
left for the retracting finger. The white curve represents the trace of the tip as defined by i2.3V . 

the course of one revolution of the spiral. To estimate this region we describe the 
boundary of the core as the trace of the tip during one revolution where the tip is 
defined as the intersection of the two contour lines 

(2.3) u* = 0.5 and v* = ~ - b , 

of the activator and the inhibitor, respectively [2J. The value v* solves J-{u* , v*) = 
with u* — 0.5. This definition ensures that the distance of the tip to the centre of 
the spiral is close to minimal. Note that this definition is arbitrary. In the large core 
limit, however, the value of the core radius r c calculated as the radius of the circle 
traced by the spiral tip does not vary much proportionally for different choices of the 
contour lines used to define the location of the tip. 

3. Freezing method. In this Section we briefly outline the so called freezing 
method as introduced in [6] before we propose our modifications to this method in 
Section |U 

We consider reaction-diffusion systems on the plane, 
(3.1) U t = DAU + f{U)=:F(U), 

with U(z,t) = {U u ---,U d ) T G K d , d > 1, z = {x,y) T G M 2 , / : R d R d and a 
diffusion matrix D G M. dxd with constant coefficients. The system (|3.ip is equivariant 
under the action of the special Euclidean group SE(2) = S 1 k R 2 consisting of ro- 
tations and translations. The Barkley model (|2.Ij) described in the previous section 
belongs to this class of equations with d = 2. 

Equivariant systems of the form (|3.1j) can be solved by the freezing method in- 
troduced in [5] for equivariant systems. This beautiful numerical method utilizes the 
fact that the dynamics of equivariant systems can be decomposed into two parts: the 
dynamics on the symmetry group and the dynamics orthogonal to it. Equivariant 
systems can be cast into a skew product form whereby the dynamics on the symme- 
try group is driven by the so called shape dynamics which is orthogonal to the group 
dynamics. This idea was developed in the seminal paper [5] and then put on a more 
rigorous footing in [UJ [TS1 [TTJ [371 [3S1 13H1 IB 131] - The method can be applied for 
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any symmetry group. Here wc restrict the exposition of the method to equivariancc 
with respect to the Euclidean group, and follow closely [SI HI]- For a more general 
framework of the method we refer to [5J [7J [301 E2] • 

3.1. General Setup. Let 7 = (0,/3) r £ S 1 x I 2 be an element of the special 
Euclidean group SE(2) consisting of the angle of rotation and translation j3 = 
(/3i,/32) T - Two group elements are linked by the operation 7W o 7^) = (0W + 

e (2) ,/3 (1) + e e <i)/3 (2) ) with 



cos 8 — sin 
sin cos 



The action a of the Euclidean group on functions U G K d can be defined by 

(3.2) [a(r/)U](z) = U(Q- e (z-0)), z G M 2 . 
The group action satisfies the properties 

(3.3) o(e) = / , a( 7 (1) o 7^) = a( 7 «)a( 7 (2) ) , 

where e is the unit element of SE{2) and I the identity matrix. Note that system 
(|3.ip is cquivariant under the action of SE(2), i.e. 



(3.4) F(a( 7 )U) = a(j)F(U) . 

The invariance of equation (|3.1j) with respect to the Euclidean group implies that it 
is possible to construct new solutions from a given solution by applying symmetry 
operations. In particular, we can rewrite a solution U (t) as 

(3.5) a( 1 (t))W(t) = U(t), 

with W(t) G R d . By formally differentiating (|3.5p with respect to t we obtain, upon 
using the equivariance condition (|3.4[) . 



(3.6) C/ t - [a 7 ( 7 )W] 7t + o(7)Wi = a(j)F(W) , 
which can be rearranged to yield 

W t = F(W) - a (7- 1 )[a 7 (7)^] 7t 

(3.7) =F(W)- S(W,i)v, 

where we used a(7) _1 = a(7 _1 ) and set v = 74. The derivative \a 1 (pj)W\ of the group 
action with respect to 7 = (0, (3) G S 1 x R 2 can be calculated as 

[a 7 ( 7 )W> = d^[W( P -e(z - /3))H 

= -WW (g„e(z - p)) ^ e ^f (* - P)vx - VW {q^{z - /?)) e _ e K v 3 ) T , 

where 17 = ©4 = u and (1/2, i^3) T = The gradient acts on vector- valued functions 
as (VW)ij = dWi/dzj. This yields the expression 

S(W, 7)1/ : = -VVK z/i£- z + g_ e (^2,^) T 

(3.8) = (yW x - - (W x ,W y ) g- e (u 2 ,v 3 f 
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Introducing new group variables (0,a) S S 1 x M 2 by setting a = Q-q/3 and denning 
parameters Hi = V\ and (a*2;M3) T = P-e^j ^3)^ allows for an elimination of the 
group variable O in (|3.8[) . Using this transformation we write (|3.8I) as 

(3.9) S{W)n ■= {yW x ~ xW y )m - W xf x 2 - W yf i 3 . 
The equations for the group variables (0, a) are given by 

(3.10) 9t = /xi with 9(0) = 0, 

(3.11) a t = Mi£f a + (M2jM3) T with a(0) = . 

For constant /ij equation ()3.11j) describes a rotation on a circle with centre at 

(3.12) {xM> y M)= (-^^ 



Ml Ml 



and radius of rotation 



(3.13) r p = J (x. p - x M ) 2 + {y P ~ Vm) 2 

for some point (x p ,y p ) on the circle. The core radius r c can be determined by apply- 
ing equation (|3.13[) to the tip of the spiral (as defined by (|2.3p ). using the centre of 
the core (xm,Vm) given by (|3.12[) . The core radius can thus be calculated from the 
group parameters, which are calculated in the freezing procedure. 

So far, the path j(t) on the group in (|3.7p is arbitrary. Therefore three additional 
degrees of freedom, equaling the dimension of the group SE(2), exist. To close the 
system we fix the location on the group orbit, and augment the equations by a phase 
condition 

(3.14) *(W,7) = 0, 

where the functional if? maps into R 3 . The phase condition vff can be chosen in several 
ways. For the time-dependent problem (|3.7p we will use the condition 

(3.15) *«in(Wi7)=/ S(W) T W t dxdy, 

JR 2 

which assigns the location on the group orbit by minimizing the temporal change of 
I [ Wt|| 2 at each time step. Note that this condition is equivalent to requiring that W t is 
orthogonal to the group orbit at W. This condition allows us to determine the freezing 
parameters fii, which are implicitly contained in (|3.15[) through Wj. These values are 
then subsequently fed into the shape dynamics to update W (see Section 13.2.11 for 
details on the implementation). 

When considering the stationary problem of (|3.7[) . we will use the following slightly 
simpler phase condition 

(3.16) *fix(W)= I S(W ) T (W -W )dxdy, 



which determines the location on the group orbit by minimizing the distance of W 
to some template function Wq. The phase condition (|3.16[) is independent of //, and 
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can therefore not be used directly in the time-dependent setting with a semi-implicit 
Crank-Nicolson scheme. Note that in order to use ^fix the template function Wo 
must be sufficiently close to a solution of (|3.1[) . For further details on phase condi- 
tions see [3 [16] . 



We summarize the closed system of partial differential algebraic equations (|3.7|) , (13. 9]) 
(|3.11|) and (|3. 14|) . The group dynamics is given by 

(3.17) e t =/ii, 9(0) =0, 

(3.18) a t = nig%a + (/i 2 ,^3) T , a(0) = 0, 
which is driven through the parameters fj,{ by the shape dynamics 

(3.19) W t = DAW + f(W) - S(W)n , W(0) = W , 

(3.20) = *(W,7), 

with S(W)fi given by ()3.9j) . This system is called the frozen system since for rigidly 
rotating spirals its solution W evolves into a stationary solution for t — > oo. 

Our aim is to numerically determine the shape of rigidly rotating spiral wave solutions 
of the Barklcy model (|2.1[) and the corresponding core radius and rotation frequency. 
Note that the shape dynamics (|3.19l) entirely determines the solution W and the 
parameters jj,i which can then be used to determine the rotation frequency oj = Hi 
and the centre of the core via (|3.12[) . To this end, it would be sufficient to solve the 
stationary problem corresponding to (|3.19[) , using for example a Newton solver (as 
done in [42]). However, the high-dimensionality of the stationary problem, needed 
for an accurate resolution of the solution, requires a sufficiently good guess for the 
initialization of the Newton-Raphson method which otherwise would not converge. 
The initial guess for the Newton solver will be generated by the application of the 
freezing method for the time-dependent problem (|3.19l) . 

3.2. Discretization. We implement the freezing method for the Barkley model 
(|2.1j) with W = (it, v) T for Cartesian coordinates and for polar coordinates. The 
choice of the coordinate system depends on the parameter range and the specific sit- 
uation. Cartesian coordinates are better suited for the investigation of the large core 
limit where a fingcr-likc solution can intersect the boundary almost perpendicularly 
and Neumann boundary conditions are a good approximation for the unbounded do- 
main. Polar coordinates are better suited for small core radii where several revolutions 
of the spiral are usually inside the computational domain. 

Simulations in Cartesian coordinates [x, y) are performed on a rectangular domain 
[0, La,] x [0, Ly] with M X N points and gridsizes Ax = L X /(M — 1) and Ay = 
L y /(N — 1). The shape dynamics reads for Cartesian coordinates as 

(3.21) d t u = Au + T(u, v) + Vu fug^ (x,y) +(/i 2 ,/i 3 ) 

(3.22) d t v = u - v + \7v fag* (x, y) T + (/i 2 , Hsf 

(3.23) = v, n) , 

with J-(u,v) = u(l—u)(u—(v+b)/a)/e. The group variables (0,a) can be determined 
by equations (|3.17[) . 
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Simulations in polar coordinates (r, if) with (x,y) = (r cos(tp), r sin(yj)) are per- 
formed on a rectangular domain [0, R] x [0, 27r — Aip] consisting of M x TV grid points 
with gridsizes Ar = R/{M — 1) and Aip = 2-k/N. The analogue to (|3.21[) reads as 



(3.25) dtv = u — v + fj,xv v + (v r , -v v )q- v (fi2,na) , 

(3.26) = *(u, «,//), 

and the dynamics of the group variables is again given by (13.171) . 

The temporal integration for both (|3.21j) and (|3.24j) . is done in discrete time steps 
At. We denote the values of the fields u and v at the n-th time-step t = nAt and 
spatial location ((i - I) Ax, (j - I) Ay) (or ((i - l)Ar, (j - l)A<p)) with i = l,...,M, 
j = 1,...,N by u™j and u"., respectively. Spatial derivatives are evaluated using 
second-order central differences. 

3.2.1. Time-dependent freezing. We use a second order semi-implicit Crank- 
Nicolson scheme to solve the time-dependent problem (|3.21[) (or (|3.24[1 ) whereby the 
linear terms are treated implicitly and the nonlinear term T(u, v) is treated explic- 
itly with an Adams-Bashforth scheme [3U [IT] . The nonlinear freezing term can be 
rendered as an effectively linear term S(— — 2 +w to be included into the Crank- 

Nicolson part of the temporal discretization, by first obtaining fi n+1 through the phase 
condition 



The additional computational costs of solving the semi-implicit equations is far out- 
weighed by the less restrictive time-step required for the semi-implicit method when 
compared to an explicit Euler method. 

In [6] an explicit Euler scheme with Neumann boundary conditions was used which 
exhibited spurious spatial oscillations in both, the interior and at the boundary of 
the domain. An upwind-downwind scheme was introduced to control the oscillations 
with partial success for large excitabilities at small values of e only. The semi-implicit 
scheme does not exhibit oscillations in the interior of the domain. However, spatial 
oscillations at the boundary persist. In Sections 14.21 and 1 4.31 we will present an expla- 
nation for the oscillations at the boundary, and provide a simple method to eliminate 
them. 

3.2.2. Stationary freezing. Travelling waves and spiral waves are both ex- 
amples of relative equilibria [35J [THj. Relative equilibria have the representation 
U(t) — a(j(t))W with a time-independent W. Hence the temporal evolution of a 
solution can be described entirely by the dynamics on the group. We may therefore 
find the spiral wave solutions of (|3.21[) by solving its associated stationary problem 



(3.24) d t u 



U rr + ~Ut + —U vv + F{u, v) + [l\U^ + (u r , -U v )g~ v (fl2, fJ^) 



(3.27) 




V 



*flx(«,«) 
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or, analogously for polar coordinates and (|3.24[) 



(A r . v u + T{u, v) + \ixu v + (u r , Iu v )q^ v (fi 2 , M3) 
u-v + mv v + (v r , Iv v )q- v {^2,^f 

with A riV u = u rr + \u r + -p2-u vv . We solve the system for the 2MN + 3 unknowns 
um,Ni vi,i, •••) Vm,N, A*i) A*2, A*3)- Typically we use M > 200 and N > 200 for 
Cartesian coordinates and M > 150 and N > 500 for polar coordinates. To solve 
this high-dimensional nonlinear system we use the Newton-Raphson method with 
line searches and backtracking (see for example |34] ) with a terminating condition 
0.5 F T F < 10~ 10 . 

4. Boundary Conditions. In the introduction, two issues were described which, 
so far, prevented a successful stable application of the freezing method to excitable 
media. Firstly, numerical instabilities in the form of spatial oscillations spoiled results 
in the case of a non-diffusive inhibitor. In the previous Section we have eliminated spa- 
tial oscillations in the interior of the domain by using a semi-implicit Crank-Nicolson 
scheme in the time-dependent freezing problem p.21|) (or p.24|l ). Spatial oscillations 
near the boundaries, however, persist. Similarly, solving the stationary freezing prob- 
lem p.28j) (or (|3.29[) ) leads to spatial oscillations near the boundary. Secondly, the 
usual Neumann boundary conditions do not reproduce the correct shape of a spiral 
in an unbounded domain. These problems are associated with the type of boundary 
condition used, and how they are formulated in the discretization. 

Naturally, computations are performed on a bounded domain and appropriate bound- 
ary conditions have to be chosen. The boundary conditions should preferably reflect 
the nature of the investigated system and its solutions. We are interested here in 
approximating spiral wave solutions in unbounded domains. By unbounded domains 
we mean either the infinite limit, where spiral wave solutions do not decay at infinity, 
or finite domains with boundary conditions but where the computational domain is 
much smaller than the actual physical domain and the physical boundaries can be 
ignored. In these cases the usual Neumann boundary conditions are, in general, not 
well suited. We discuss their impact on the freezing method in Scction [4.1l In Section 
14.21 we present spiral boundary conditions and show how they can be implemented 
for freezing methods. The standard implementation of both these boundary condi- 
tions leads to oscillations of the inhibitor near the boundary. In Section 14.31 we will 
therefore introduce an implementation for both Neumann and spiral wave boundary 
conditions, which does not exhibit spurious spatial oscillations. 

4.1. Neumann boundary conditions. Most work on spiral waves in excitable 
media uses either Dirichlct or Neumann boundary conditions (NBC) (e.g. [2]). These 
boundary conditions are physically meaningful for simulations of excitable media on 
bounded domains, e.g. in chemical experiments. However, if simulations are per- 
formed with the intention to understand the behaviour of spirals in unbounded do- 
mains, they have the disadvantage of not respecting the underlying symmetry at the 
boundary. Nevertheless, Neumann boundary conditions have been used extensively. 
For simplicity we restrict our discussion to polar coordinates, for which Neumann 
boundary conditions are formulated as u r = and v r = at r = R. These are 
discrctized according to 

(4.1) Um+1,3 = "M-1,3 and v M +i,j = «M-i,j , 
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for j = 1, • ■ • , N. The values of um+ij, VM+i,j can then be used in the evaluation of 
the diffusion and advection terms on the boundary. 




(g) SBC for u, and (PS) (h) (i) 

for v 



Fig. 4.1. Inhibitor v of a frozen spiral wave solution calculated via the stationary problem 
\S. 2.91 ). From left to right: representation in the Cartesian plane (x,y) = (r cos((p),r sin(tp)), in the 
polar (r,ip) -plane, and a close-up of the radial boundary, (a)-(c): Neumann boundary condition 
(NBC) {4- ljj , (d)-(f): spiral boundary condition \4- f or both, activator and inhibitor, (g)-(i): 
spiral boundary condition (SBC) {4. 76 for activator only and free boundary condition \4-23ty for 
the inhibitor. We chose e = 0.025 in the Barkley model 12. 16 . and R = 21.74, Ar = 0.1257 and 
Aip = 0.01 for the spatial discretization. 



In the two top left panels of Figure 14.11 we show contour plots of the inhibitor 
v calculated by solving the stationary problem (|3 . 29[) in polar coordinates. Figure 



4.1(a) shows the spiral solution on the circular domain with radius R in the Cartesian 



plane (x,y) = (r cos(y), r sin(<p)). In Figure 4.1(b) we show the same solution in the 
(r, </?)-plane. One sees clearly how the contours bend towards the boundary to satisfy 
the Neumann boundary condition u r = 0, v r = at r = R. This kink is localized 
near the boundary and does not extend far into the domain. However, the bending 
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of the spiral wave solution leads to an effective wider transversal cross-section of the 
solution near the boundary. A wider activator profile allows the inhibitor to adopt 
larger amplitudes (cf. the darker color (online red) of the inhibitor in the close-up in 
Figure 4.1(c) I. The boundaries, however, do not affect the dynamics near the core - 
provided they are located sufficiently far away from the spiral wave tip - and hence 
the presence of a kink in the solution near the boundary does not affect the values of 
the rotation frequency ui and the core radius r c . 

More important, from a numerical stability point of view, is the following issue. 
In the close-up of the inhibitor near the boundary in 4.1(c) strong spatial oscillations 
are clearly visible. For a more detailed view we show a slice of the inhibitor near the 
boundary along the grid line ip — 3tt/A in Figure 4.2(a) These oscillations are a well 
known problem of the freezing method for excitable media, and are associated with 
the non-diffusive nature of the inhibitor [6j |4TJ [8j . These spatial oscillations occur 
only in the non-diffusive inhibitor v, and are absent for the diffusive activator u. We 
have checked that there are no oscillations for sufficiently large diffusion coefficient 
D v 7^ of the inhibitor. 

The amplitude and the extent of the oscillations increase strongly for larger values 
of e, which prohibits an accurate numerical analysis of the large core limit. 




2/3 r / R 1 2/3 r / R 1 

(a) (b) 

Fig. 4.2. Left: Slice of the close-up in Figure \j. 1 (c)\ along a radial grid line for the inhibitor 
solution v where NBCs are used for activator u and inhibitor v, demonstrating spurious spatial 
oscillations near the boundary. Right: The same slice but now with SBC for the activator u and the 
one-sided derivative \4-23j) as boundary condition for the inhibitor v, corresponding to Figure\4-l(i)\ 



4.2. Spiral wave boundary conditions. In this Section we introduce two 
boundary conditions based on simple geometrical approximations of spiral waves in 
order to mimic the behaviour of spirals in unbounded domains, or in finite domains 
when the size of the computational domain is much smaller than the physical domain. 
A natural version of this type of boundary conditions using Archimedean spirals was 
introduced in jlOj . Therein, spiral wave solutions were "grown" with a predefined 
wavelength A. We will expand and generalize this idea with the aim to apply it to 
the freezing method. Here the wavelength is not known a priori and, moreover, the 
centre of the spiral wave generally moves during the freezing procedure. In addition 
to the approximation by an Archimedean spiral, we will present boundary conditions 
where the spiral wave is approximated by an involute of a circle [45] . The two ap- 
proximations coincide in the far field but differ near the spiral wave tip. We will 
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formulate spiral boundary conditions for polar and Cartesian coordinates. It has to 
be noted that both, the Archimedean spiral and the involute of a circle, are just 
geometrical approximations of contour lines of the actual spiral wave solution in un- 
bounded domains. There exists, to our knowledge, no rigorous theoretical justification 
for these approximations. However, we will see that these approximations may serve 
as convenient constructs to formulate boundary conditions. 

We assume a spiral wave solution U(f, ip) with a constant far field wavelength A, 
centred at the origin (xo,yo) 01 the polar coordinate system (f , (p). A common choice 
to describe such a spiral wave geometrically is by means of an Archimedean spiral. In 
this approximation contour lines of the spiral wave are given by 

(4.2) $ s = fhf + (p = const , 

allowing us to simplify U(f, tp) = V s (fhf+ip). The parameter fh is given by fh = 27r/A, 
which can easily be seen by noting that a spiral wave profile is 27r-periodic with 
V s (<b s ) = V s (<t s + 2ir), and therefore fhf + (p = fh(f + 5^) + (p. The approxima- 
tion of a spiral wave solution of an excitable medium by an Archimedean spiral is 
fairly accurate in the far field, away from the spiral wave core, but fails close to it. 
This is illustrated in Figure l4~3[ where two examples of rigidly rotating spiral wave 
solutions of the Barkley model (|2.1[) are shown for different excitabilities e. The spi- 



ral wave in Figure 4.3(a) rotates around a circular core which is small compared to 
the computational domain. The contour lines of this spiral wave coincide well with 
a fitted Archimedean spiral (light dashed line; online: green). In comparison, for 
higher values of e, when the core of the spiral wave solution is larger, as depicted in 



Figure 4.3(b) (smaller dashed circle; online: red), the contour lines of a spiral wave 
solution are approximated by an Archimedean spiral (light dashed line; online: green) 
only sufficiently far away from the core. This effect worsens for larger core radii. The 
inability of Archimedean spirals to approximate spiral wave solutions near the core 
can be understood by considering that spiral waves possess locally an approximately 
constant velocity which is normal to the wave front. However, for uniformly rotating 
Archimedean spirals the radial velocity is constant. This is particularly problematic 
close to the origin where the normal direction of the wave is considerably different to 
the radial direction. Especially for numerical investigations in the large core limit, 
where numerical domains are not able to include several spiral wave coils, Archimedean 
spirals will not serve as good approximations within the computational domain. 

Note that logarithmic corrections to (|4.2j) have been considered (see for example 
[55]). However, these corrections become negligible in the far field near the boundary, 
and more importantly for our purpose here as boundary conditions, their derivatives 
with respect to the f will be small at the boundary. 

Based on the assumption of a locally normal velocity for spiral waves, it was suggested 
in [HI H5J [5S] to approximate spiral waves by involutes of a circle. In this case contour 
lines of the spiral wave are given by 



(4.3) = p + arccos(rj/f) — \J (f/ri) 2 — 1 — — = const , 



allowing us to simplify U(f,<p) = Vi(cp + arccos(r//f) — \J (f/ri) 2 — 1 — 7r/2). The 
parameter 77 denotes the radius of the circle that is revolved by the involute. The 
involute is only defined for f > 77. The radius 77 is somewhat arbitrary and we have 
to choose a proper definition that suits our application. We will later explain how to 
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(a) e = 0.025 (b) e = 0.065 

Fig. 4.3. Spiral wave solution u calculated by solving the stationary frozen system H3.29V on 
a circular domain of radius R = 21.74. (a) Spiral solution for e = 0.025 with a superimposed 
Archimedean spiral (dashed line; online: green) with fh = 2ir/\ = 0.59. (b) Spiral solution for 
e = 0.065 with superimposed Archimedean spiral (light dashed line; online: green) with in = 2tt/X = 
0.13 and superimposed involute of a circle (dashed black line) with rj = 7.9 (dashed circle; online: 
magenta). The smaller circle (solid line; online: white) with radius r c = 5.3 is the trace of the tip 
of the spiral wave u defined in 12.36 . 



choose an appropriate rj. Note that for r — >• oo we have $ s — — > 0, i.e. the two 
approximations of Archimedean spiral and involute of a circle coincide in the far field. 



In Figure 4.3(b) we show how the approximation of an involute of a circle compares 
to the one by an Archimedean spiral. In |33| the two geometric approximations 
were compared in their ability to approximate experimental data of the Belousov- 
Zhabotinsky reaction. As the experimentally obtained spiral waves were rotating 
around small cores, both approximations showed equally good agreement. It is clear 
that both approximations fail close to the spiral tip |50j . however, for our purpose, 
as illustrated in Figure 4.3(b)| involutes of circles are better suited, especially for the 



large core limit. 

We can use the approximations of an Archimedean spiral or of an involute of a 
circle, to formulate boundary conditions. Under these approximations we can use 
contour lines to express spiral wave solutions as 

(4.4) U(r,(p) = V s , I ($ s , I (r,®), 

where $ Sj j(f, <p) is given by equation (|4.2[) or (|4.3p . respectively. Differentiating (|4.4I) 
leads to 

(4.5) U f = &(r, 0) U$ , 

where a(f, ip) = <9$ s ,/(f, <p)/dr is a coefficient depending on which spiral approxima- 
tion has been chosen. For the Archimedean spiral we find 

(4.6) a(f, (p) = m . 

Note that <5 is constant for Archimedean spirals. For the involute of a circle we find 



(4.7) = ~\l 
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When evaluated at the boundary of the computational domain, we coin this type of 
boundary condition spiral boundary condition (SBC). SBCs have the advantage that 
derivatives of a spiral wave solution U on the boundary can be expressed entirely by 
known values of U from inside the domain and from the boundary. Note that the 
geometric approximations can still be used to formulate SBCs, even, if close to the 
tip, they are not accurate. 



The spiral boundary conditions (|4.5j) are formulated within the coordinate system 
(r, tp) which is the coordinate system with origin at the centre of the spiral wave 
{xo,yo)- This does, in general, not coincide with the polar grid of the computational 
domain. Assume computations are performed on a circular domain of radius r = R 
in a coordinate system (r, tp) with centre (0, 0). 

To complicate things, the centre of the spiral {xo,yo) arj d the associated coordi- 
nate system (f, ip) typically shift during the process of freezing. We therefore need to 
perform a coordinate transformation relating the two coordinate systems, (r, tp) and 
(r, tp) , in order to express the spiral boundary conditions ([4.5)1 in terms of the coor- 
dinate system of the computational domain (r,tp). Using elementary trigonometric 
relations (see Figure 14.4)1 we can write 



(4.8) r = \l (r cos tp - x ) 2 + (r sin tp - y ) 2 , 



(4.9) tp = arctan 



r am tp- yo 
r cos tp — Xq 



We can now formulate SBCs on the actual computational domain by inserting the 
transformation into (14. 511. to obtain 



(4.10) U r = a{r,tp)U v , 
with 

(1 a(r,tp>)r cos(p - ip) + sin (93 - tp) 

(4.11) ot{r,tp) = - 



r —dt(r, tp) f sm(tp — tp) + cos(p — tp) 



where <S(r, tp) is given by ([4.6)1 for Archimedean spirals and by (|4.7[) for involutes of 
circles. 



4.2.1. Determination of the parameters of the spiral boundary condi- 
tion. The SBC requires knowledge of the location of the centre of the spiral wave 
(^0) Vo) and its wavelength in the case of the Archimedean spiral (|4.6|) . and knowledge 
of the parameter 77 for the case of the involute of a circle (|4.7[) . In the following we 
will explain how these parameters can be determined. 

We start with the Archimedean spiral. The freezing method automatically deter- 
mines the instantaneous centre of rotation of the spiral wave solution via p,12[) . We 
may therefore determine (xo,yo) on the fly during the freezing procedure as 

(4-12) (xo,yo) « (xm,Vm) =( -,— ) , 

V Mi Mi/ 

which becomes exact once the group parameters p,i have become constant. Note 
that for rigidly rotating spirals /jj = u ^ unless we are at the bifurcation from 
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Fig. 4.4. Coordinates of a point (x p ,y p ) in two different polar coordinate systems (f , (p) and 



rigidly rotating spirals to retracting fingers. Approaching this limiting case, (im, Vm) 
becomes larger, reflecting the increase of the core radius which diverges at criticality. 
If the tip were fixed at the origin, the core radius could be determined exactly by 
r c = \\(xm, Um)\\ using only the group parameters. In this case the rotation frequency 
and the core radius are exactly inversely proportional to each other. In the large 
core limit, an offset of the spiral tip from the origin of the computational domain, 
however, becomes negligible as demonstrated later in Figure 5.7(b) See [TB] where a 
phase condition is used which pins the tip at the centre of the domain. 

The parameter rh = 2n/X can be estimated on the fly as well during the freezing 
procedure. For small core spirals we can do so by employing the dispersion relation of 
travelling wave trains [55]. Before the application of the full 2D-frcezing procedure, 
the dispersion relation c w (X) for the velocity of a travelling wave train is determined 
for each fixed e. This is done by applying the freezing method to the Barkley model 
(12. f p in a ID periodic domain with length A. Assuming that the ID velocity c w (X) 
is a good approximation for the normal velocity of the far field spiral wave coils, this 
velocity should equal the asymptotic far field velocity of a spiral wave c s = (ui/2tt)X, 
which is consistent with our geometric approximation of Archimedean spirals. The 
rotation frequency u> is determined during the freezing procedure and given simply 
by u; = The wavelength A of the spiral wave is then obtained as the solution 
of c s (A) = Cm (A). For sufficiently large core radii, when the spiral wave coils do not 
interact with each other, we may replace c w (X) by the ID velocity of an isolated pulse 
Coo. The wavelength can then be estimated by A = 2-KCoajuj. The velocity Coo can 
also be determined using a ID freezing procedure, or via direct simulations of the ID 
version of the Barkley model (|2.1j) . where the box length is chosen large compared to 
the decay length of the inhibitor. 

For large core spirals, a simpler method can be used to determine A. In this limit 
we can employ the approximation A = 2irr c which is independent of Coo or the dis- 
persion relation c w (X), which we would need to determine in advance. In Section l5.II 
we will present numerical results corroborating these approximations. 
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The SBC involving the approximation of involutes of a circle requires that we specify 
the radius 77 of the circle around which the tip revolves. We found that the involute 
of a circle with radius 77 = c^/lo is a better approximation of contour lines of the 
actual spiral wave solution than the involute of a circle with the somewhat arbitrary 
radius r c , in particular in weakly excitable media with large core radii. 

4.2.2. Applicability of the spiral boundary conditions. In the following 
we investigate the applicability of the spiral boundary conditions (|4. 10|) and (|4.1ip . 
There are two possible problems that can occur: Firstly, the case when there are grid 
points on the boundary at which U v = but U r ^ 0, violating (|4.10[) . Secondly, the 
expression for a(r, <p) in (|4.11|) can possess singularities. We will derive conditions on 
the location of the spiral centre (xq, yo) to ensure that the spiral boundary condition 
is applicable and both possible problems are avoided. For simplicity, we will use the 
example of the involute of a circle. 

At first we investigate, when there are points on the boundary such that U v = 
but U r 0. This is equivalent to asking whether there exist points on the boundary at 
which the involute is tangent to the circular boundary of the computational domain. 
The involute of a circle with radius 77 can be described by 

(4.13) I(s) = 77 ( cos ) s + s sins sin s — s cos s + ( x ) yo , s <G [0, 00) . 
A parametric equation for the circular boundary reads as 

(4.14) C{a) = i? ( cos )ct shier, cre[0,27r). 

Without loss of generality we may set yo — 0. A necessary and sufficient condition 
for a tangency at the boundary is given by the solution of the two equations C T> t = 
and C = I, which is found to be xq cos s = — 77. Hence, SBCs are violated whenever 

No I > Ti. 

In a next step we investigate the conditions for singularities of a(r,tp). Zeros of 
the denominator can be found for tan(^j — (p) = \/(r a{f, ip)). Using the trigonometric 
identity tan(<£ — tp) = (t&ntp — tany>)/(l + tantptancp) and the coordinate transfor- 
mation (|4.8I) . (|4.9[) we find xq simp = 77, implying as before |xo| > 77. We note here 
without derivation, that for Archimedean spirals one finds |a;o| > r c [24] . 

In the small core limit, when the core radius is small compared to the computational 
domain, the centre of the core can be placed close to the origin of the computational 
domain, assuring \xo\ < r c . If the core becomes larger upon increasing e, the centre of 
the spiral wave core will move outside the computational domain, if the spiral wave 
tip remains resolved within the domain. In this case the conditions \xq\ < rj can 
be satisfied, by placing the tip of the spiral wave solution close to the centre of the 
computational domain, which implies |xo| ~ r c . Recalling our definition 77 = Cqo/u;, 
we have 

(4.15) r/ = £2£>££ =rc; 

Id UJ 

where we define c c to be the velocity of the spiral tip tangential to the circle with 
radius r c . Because of the positive curvature of the wavefront close to the core, this 
velocity is smaller than the velocity Coo of the spiral coils in the far field with vanishing 
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curvature. Hence, if the tip is placed close to the centre of the computational domain, 
\x \ < rj is satisfied, assuring that the SBC is neither singular nor violated. 

In the large core limit, however, several problems arise, which cannot be resolved 
by an appropriate positioning of the spiral wave solutions within the computational 
domain. For the approximation with involutes of circles different problems may arise 
which need to be addressed. First, for large core spirals the involute circle (or the 
spiral core) intersects the boundary of the computational domain, and therefore there 
are parts of the boundary for which f < rj, and the SBC (|4.10[) with (|4.11j) and 
(|4.7I) is not defined anymore. However, in this case we can formulate mixed boundary 
conditions. We keep SBCs on the part of the boundary where f > 77, and we impose 
the NBC U r = for the remainder of the boundary where r < fj, i.e. 



(4.16) U r 



a(r, ip)U v for f > n 



for f < n 



where a is given by (|4.11[) using (|4.6p for the Archimedean spiral or (|4.7I) for the 
involute of a circle. This is justified since, per definition, the core region of a rigidly 
rotating spiral is the part of the domain which is never excited by the spiral wave, 
with U ~ and U r ~ 0- For the part of the domain which does not lie within the 
core region, but within the circle of radius 77 > r c , we may also set U = U r = since 
in the large core limit, the fields decay fast enough. 

Approaching the critical point with r c —> 00, we inevitably reach the point when 
77 — r c is larger than the domain size R. At this point none of the geometric approx- 
imations we discussed is valid anymore. In this case, NBCs, formulated in Cartesian 
coordinates, prove to be a good approximation, as the curvature of the spiral wave 
solution becomes negligible and the spiral wave appears to have the shape of a trav- 
elling finger. 

Analogously to (|4.10j) we can formulate spiral boundary conditions for Cartesian grids 
(x,y). We perform a transformation between (f, (^-coordinates of the spiral system 
and (x, y)-coordinates of the computational grid 



(4.17) r = (x ~ xof + (y ~ yof , 

(4.18) (p = arctan ( y ~ y ° ) . 

\X-XqJ 

We find 

(4.19) U x =a{x,y)U v 
for boundaries parallel to the y-axis, and 

(4.20) U y = a(x, y y 1 U x 

for boundaries parallel to the £-axis. Here a(x,y) is given by 

a(f, ip) f cos if — sin (p^ 



(4.21) &{x,y)- 

air, tp) r sin ip + cos tp 

where a is given again by (|4.6j) for the Archimedean spiral and by (|4.7[) for the 
involute of a circle. The variables f and Cp are expressed as functions of (x, y) using 
the coordinate transformation (|4.17[) . 
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Spiral boundary conditions in Cartesian coordinates can only be applied for large 
core spirals when the core radius is considerably larger than the box length of the 
rectangular computational domain. This is due to the fact that on rectangularly 
shaped domains covering at least one revolution of the spiral arm, there are always 
points on the boundary at which the spiral arm is tangent to the boundary, violating 
the SBCs (|4.19|) or (|4.20|) . Therefore we have to restrict the method to large core 
spirals where only a finger-like part of the spiral is resolved within the computational 
domain. The orientation of the domain can be chosen such that only one boundary 
intersects with the spiral wave arm. In this case we can set-up boundary conditions, 
by imposing an SBC on that boundary and Neumann boundary conditions on the 
remaining boundaries where activator and inhibitor are assumed to be approximately 
zero. Without loss of generality we choose y = to be the boundary which intersects 
with the spiral. The mixed SBC-NBC boundaries are then written as 

for x = f a(x, y)~ x U x for y = 

and U y = I 

for x = L x I for y = L y 



In Figure 4.5(a) a contour plot of the activator u in Cartesian coordinates is shown, 



where we used boundary conditions (|4.22|) . Figure 4.5(b) shows an overlay of this 
solution with a similar solution where computations were performed with NBCs for 
all boundaries. We show the contour lines of u = 0.5 and v = 0.5 a — b respectively, 
and an inset zooming into the area close to the boundary. For the solution obtained 
by using SBCs, the boundary appears transparent, whereas the solution obtained by 
using NBCs for all boundaries exhibits a kink near the boundary. The difference 
between the two solutions is confined near the boundary, and is absent in the interior. 





25 50 15 20 25 

X X 

(a) (b) 

Fig. 4.5. Spiral wave solution u in the large core limit obtained by solving the stationary frozen 
system 13.286 in Cartesian coordinates. The excitability parameter is e = 0.0799, and the spatial 
discretization is Ax = Ay = 0.125. (a) Contour plot of activator u. The inset is added to show the 
correct aspect ratio, (b) Contour lines u = 0.5 (online: black, green) and v = 0.5a — b (online: blue, 
red) for SBC (involute of a circle) and NBC respectively (inset: zoom into dashed box close to the 
boundary at y = 0). 



4.3. Oscillation- free implementation of boundary conditions. In the mid- 
dle row of Figure I4TT1 we show results of the stationary problem (|3.29[) in polar coordi- 
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nates using the spiral wave boundary condition (|4.10p for the involute of a circle with 



(|4.11[) and (|4.7I) . Comparing Figure 4.1(b) and Figure 4.1(e) we see that the spiral 



boundary condition respects the shape of the solution at the boundary and does not 
include spurious kinks. However, as can be seen in Figure 4. 1 (f ) [ the spiral boundary 



conditions as described in this section, are not able to control the oscillations near 
the boundary. 

In this section we present an explanation for these oscillations, and suggest a 
simple way to eliminate them. To understand what causes spurious oscillations for 
NBC and SBC, we investigate their numerical implementation in more detail. Without 
loss of generality we restrict the discussion here to polar coordinates. 

For both boundary conditions, NBCs and SBCs, the equation for the inhibitor 
involves only first derivatives. This implies, that at the radial boundary (R = 
(M — l)Ar, (j — 1)A(/?), the inhibitor v is computed only from values of v and u on the 
boundary from the previous time step. Whereas the discrete Laplacian, present in 
the equation for the activator u, couples values Umj at the boundary to those in the 
interior, i.e. um—xJj the boundary values of the inhibitor Vmj arc decoupled from 
the interior, and do not receive the outward flowing information from the interior. 
It is this decoupling of the boundary from the interior for the non-diffusive inhibitor 
which causes the spatial oscillations near the boundary. (We recall that there are 
no oscillations for the diffusive activator u, and also no oscillations when diffusion is 
added to the equation for the inhibitor.) 

For both cases, NBCs and SBCs, we propose a simple method to overcome this de- 
coupling. For the activator u we use NBCs or SBCs as discussed. However, instead of 
invoking this boundary condition for the inhibitor as well, we evaluate the derivative 
of the inhibitor at the boundary by a one-sided second-order discretization according 
to 

(A on „ _ 3v MJ ~ 4VM-1,3 + VM-2,j 

^■ Z6 > v r\M, } - 2Ar 

For Cartesian coordinates we use equivalent expressions. 

Considering the original unfrozen Barkley system (|2.ip , there are two reasons why 
this simple trick works. Firstly, the boundary condition is consistent with the fact 
that information flows outwards when studying spiral waves, and therefore outer grid 
points are influenced by inner grid points only. Secondly, the inhibitor is "slaved" to 
the activator in the sense that if the activator is known, the linear equation for the 
inhibitor can be integrated explicitly. The boundary conditions of the activator are 
therefore (approximately) inherited by the inhibitor. 

In Figure |4~T1 (bottom row) we show how the implementation of the spiral wave bound- 
ary condition for the activator and the one-sided derivative (|4.23p for the inhibitor in 
the stationary problem (|3.29p in polar coordinates suppresses the spatial oscillations 



and avoids changes in shape and amplitude (see also Figure 4.2(b)). In Figure 
we show how the implementation of Neumann boundary conditions for the activator 
and the one-sided derivative (|4.23|) for the inhibitor controls the spatial oscillations 
and the spiral wave shape in the stationary problem f|3 . 28[) in Cartesian coordinates 
in the weakly excitable regime. 



During extensive numerical tests |24j we found that in highly excitable media only spi- 
ral boundary conditions (|4.10p in conjunction with the one-sided boundary condition 
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(|4.23[) give accurate results for the stationary and the time-dependent freezing method 
without spurious oscillations and without unnatural deformations of the spiral shape. 
For the large core limit, in weakly excitable media, Neumann boundary conditions 
in Cartesian coordinates are the optimal choice to avoid numerical oscillations and 
to properly resolve the shape of the spiral wave solution. Note that in the large core 
limit NBCs resemble SBCs. This is due to the fact that the overall curvature of the 
spiral wave, captured in the computational domain, is close to zero. In a way NBCs 
can be viewed as a simpler and more convenient formulation of SBCs in the large core 
limit, which do not involve the core radius and the wavelength of the spiral wave. 

In the following we will be using SBCs in polar coordinates for highly excitable 
media, and NBCs in Cartesian coordinates in the large core limit. We will also em- 
ploy the one-sided boundary condition (|4.23l) in all numerical simulations to avoid 
oscillations. 



0.6 




(c) (d) 



Fig. 4.6. Inhibitor v in the weakly excitable case with e = 0.075, calculated from the stationary 
frozen system (3.28\) in Cartesian coordinates with NBCs for activator and inhibitor (top panels) 
and NBCs for activator only and one-sided derivative \4 . 23]) for the inhibitor (bottom panels), (b) 
and (d) are close-ups of the respective solutions depicted in (a) and (c), zooming in on the boundary 
at y = 0. 
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5. Numerical investigation of the large core limit of spiral waves. In 

this Section we provide numerical results on the wavelength A, the core radius r c and 
the rotation frequency ui of spiral waves in the large core limit. The rotation frequency 
u> and the core radius r c are computed from the group parameters /ij. We recall that 
w = fii and r c can be calculated using (|3.13[) with the centre of the spiral given by 
(1 3 . 1 2 [) and the location of the tip defined as the intersection of the contour lines (|2.3I) . 
We approach numerically criticality where the core radius develops a singularity and 
the frequency approaches zero. Before we present these results, we collect and discuss 
some practical aspects of the different methods we introduced to study spiral waves 
in the large core limit. 



We solve the stationary freezing problem (|3.28[) . increasing e. As initial guess for 
the Newton-Raphson method we use the frozen solutions of the previous value for e. 
We use a discretization of Ax = Ay — 0.125 for Cartesian simulations, unless stated 
otherwise. We have checked the quadratic convergence of the error in the calculation 
of the core radius and the rotation frequency with respect to the spatial discretization, 
and found that at this resolution the values have sufficiently converged. 

For polar coordinates we use a discretization Ar = 0.125 and Aip = 27r/640 « 
0.01. For a computational domain with R = 20 this corresponds roughly to an equiva- 
lent discretization of Ax = 0.125 and Ay 0.2 near the boundary. This immediately 
alludes to practical limitations of polar coordinates for large computational domains; 
wave profiles away from the centre of the computational domain are not properly 
resolved, and R has to be sufficiently small, unless one employs a computationally 
expensive fine angular discretization. At this point it is important to repeat, that in 
highly excitable media with small core radii at sufficiently small values of e, polar co- 
ordinates work very well. This is due to two effects: First, transverse wave profiles are 
wider for smaller values of e, and second, the smaller wavelength in this case causes 
a cross section of the spiral wave at the boundary which is wider than its transverse 
profile. 



In Figure 15.11 we show results for solving the stationary frozen system (|3.29|) in polar 
coordinates. We used two types of boundary conditions, i.e. NBCs (red o) and SBCs 
based on the approximation by involutes (blue •). We further included results from 
direct simulations of the full Barklcy model (|2.1[) (black x). For a computational 
domain of size R = 21.74 direct simulations are limited to excitabilities of e < 0.065 
with corresponding radii r c < 5.3. One can, in principle, determine the core radius for 
parameter values larger than e = 0.065, however, only with the additional computa- 
tional cost of increasing the computational domain. There will be a value of e where 
this is not possible anymore with given computational power. We stress that this 
value of e is way below what we call the large core limit. This can be seen in Figure 
15.21 (left panel) where we show contour plots of the activator at different values of e 
corresponding to data points in Figure [5T] for SBCs. The white circular lines represent 
the trace of the tip as defined in (|2.3p and are calculated by solving equations ([3. 171) 
for the group variables. Figure 5.2(c) indicates that it becomes inherently difficult to 
calculate spirals and their characteristic parameters r c and u for large cores by direct 
simulations of the full Barkley system (|2.1j) . It is this observation which makes the 
freezing method so attractive for large core spirals. 

The insets in Figure I5TT1 clearly show that SBCs outperform NBCs in the highly 
excitable regime, and allow for the determination of frozen solutions and the values 
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of the core radius r c and the rotation frequency oj for a greater range of parameter 
values. The breakdown of NBCs is linked to the finite size of the computational do- 
main as illustrated in Figure 15.31 We can see that the finger is not properly oriented 
in the finite domain and that therefore NBCs are not a suitable choice, bending the 
solution into an unnatural direction. This does not happen for the shape-preserving 
SBCs. The negative effect of the (inaccurate) NBCs to deform contour lines of the 
solution is proportionally larger for moderate computational domains than for larger 
ones (which makes polar coordinates more sensitive to these effects then Cartesian co- 
ordinates). In principle, one may extend the range of validity in e-parameter space for 
NBCs by considering larger and larger computational domains; however, this would 
quickly become computationally unfeasible. 
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Fig. 5.1. Core radius r c and rotation frequency ui of spiral waves as functions of e computed 
in polar coordinates. Here 'x ' denotes values obtained by a direct simulation of the Barkley model 
12. and '•' and 'o' represent results from solving the stationary frozen system (3.29\) with NBCs 
and SBCs respectively. A computational domain with R = 21.74 was used. 



To approach the large core limit, we employ a Cartesian coordinate system for values 
of e > 0.068. In Figure 15.41 we show the rotation frequency and the core radius 
of simulations where we approach the critical point at e c « 0.08054091. Here we 
use L x = 50 and L y = 62.5. In Figure 15.21 (right panel) we show contour plots of 
the activator solution of the frozen stationary system (|3.28[) at different values of e 
corresponding to data points in Figurc [5~T1 for NBCs. The exact value of e c depends on 
the discretization and also on the size of the computational domain. See Section 15.21 
for a discussion on this issue. Compare the range of r c attainable in polar coordinates 
and Cartesian coordinates (cf . Figure 15.11 and inset of Figure I5.4|) . Contrary to the 
results with polar coordinates, in Cartesian coordinates NBCs are applicable for a 
larger range of parameter values than SBCs. As discussed in Section 14. 2. 2[ SBCs are 
not applicable in the large core limit, when the spiral wave appears as a finger in 
computational domains of finite size. The part of the spiral wave solution, starting 
at the tip, which cannot be approximated by an Archimedean spiral or an involute 
of a circle, increases the closer one is to criticality. As a proxy for the extent of this 
region we depict in Figure [5~5l how n — r c grows with e — > e c . For values e > 0.0797, 
corresponding to core radii r c > 2523, we have ri — r c > L y , and the freezing method 
using SBCs is not applicable anymore. NBCs, on the other hand, become a better 
approximation the closer we are to the critical point, where the spiral wave solution 
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Fig. 5.2. Contour plots of the activator solution of the stationary frozen system in polar 
coordinates with SBCs based on approximations by involutes (left) and in Cartesian coordinates 
with NBCs (right) for increasing values of excitability e. The white circular lines indicate the trace 
of the tip. 
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(a) SBC 



(b) NBC 



Fig. 5.3. Contour plot of the activator u found as the solution of the stationary frozen system 
at e = 0.079 with R = 21.74 with SBCs (left) and NBCs (right), respectively. 



approaches zero curvature - provided that the finger is appropriately oriented within 
the computational domain to assure a perpendicular intersection with the boundary. 

To study the behaviour of the wavelength, the rotation frequency and the core 
radius of spiral waves in the large core limit, we therefore use from now on Cartesian 
coordinates and Neumann boundary conditions. In Figure lOH we show contour plots 
of the activator and the inhibitor close to criticality, illustrating the appropriateness 
of Neumann boundary conditions in Cartesian coordinates in this limit. 
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Fig. 5.4. Core radius r c and rotation frequency ui of spiral waves as functions of e with e < e c 
computed in Cartesian coordinates. Here 'x ' denotes values obtained by direct simulation of the 
Barkley model fli.il) . and ' and 'o ' represent results from solving the stationary frozen system 
(3.28\) with NBCs and SBCs respectively. A computational domain with L x = 50 and L y = 62.5 
was used. 



5.1. Wavelength. The application of the spiral boundary conditions using Archi- 
medean spirals requires the knowledge of the wavelength A. In the small core limit 
the wavelength can be determined as the solution of the implicit equation 

(5.1) UA) = ^A, 
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Fig. 5.5. Core radius r c and radius of the circle of the involute rj = Co^/lo as a function of e. 
The critical value of the excitability parameter is e c = 0.08054091. 
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Fig. 5.6. Contour plot of a spiral wave solution of system f3. 28V at e = 0.08052, close to the 
critical excitability e c . The (a) activator and (b) inhibitor are shown with the inset depicting them 
with the correct aspect ratio. A computational domain with L x = 50 and L y = 140 was used. 



where c w (X) is the velocity of a ID wave train with wavelength A (sec Section |4. 2. 
In the large core limit, where the inhibitor decays sufficiently quickly and spiral wave 
coils do not interact, we may further simplify to A = 2-71X00/0;, where Coo is the 
velocity of an isolated ID pulse. In Figure 5.7(a) we show a comparison of these 
expressions with numerical results from a direct simulation of the full Barkley model 
(|2.f I) . The velocities c w (X) and Coo are determined by freezing pulses in the corre- 
sponding ID- model with box length A using the same discretization Ax = 0.125 as 
in two dimensions. We see that the wavelength determined by (|5.1|) is a reasonably 
good approximation of the true wavelength even for small core radii. For larger radii 
the two methods to determine A converge. 

In the large core limit we can deduce a simpler approximation for the wavelength 
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which does not require the independent determination of the ID velocities. One 
can define two approximate temporal periods for rigidly rotating spiral waves with 
wavelength A and curvature k. First, the temporal period T p = A/c(A, k) of a spi- 
ral wave with velocity given to first approximation by c(A, k), and second, the time 
T r = 2n(r c + 5)/c n (X, k) which measures the time of one revolution of a spiral wave 
tip around a circle with radius r c + S chosen such that the normal velocity c n (A, k) of 
the spiral tip is tangential to that circle. Equating these two temporal periods leads 
to the kinematic relation ED 



(5.2) 



A = 2ir{r c + S) 



c(A, k) 

Cn(X, K) 



In the large core limit, we expect c ra (A, n) = c(A, k) = c 00l and 6 <^ r c (see Figure IS~5 
In this case ()5.2j) reduces to the simple relationship 



(5.3) 



A = 2irr r 



In Figure 5.7(b) we show that (|5.3[) is a good approximation of the wavelength in 
the large core limit and matches well with A = 2ttCoo/uj. Figure 5.7(b) shows that 



U) ~ r' 1 as already indicated by (|4.12[) implying n\ + n\ 



3 , resembling the motion 



of travelling waves. 
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Fig. 5.7. Different estimates of the wavelength X. Red 'o': A obtained by using the nonlinear 
dispersion relation of ID wave trains c w (X) = (ui/2tt)X. Blue X = 2-kCco/oj. (a) Small core 
limit. Black 'x ': Values obtained from direct simulations of the full Barkley model 12. (b) Large 
core limit. The dashed reference line corresponds to X = 2nr c . Here ui and r c are obtained by solving 
the stationary frozen system [3. 286 . 



5.2. Scaling behaviour of the large core limit. In Figure lS~8l we show results 
on the scaling behaviour of the rotation frequency tu and the core radius r c as a 
function of the distance to criticality (e — e c ). We can clearly identify a linear scaling 
regime 

(5.4) lo - (e - ec) 1 and r c ~ (e - e c ) _1 

at the bifurcation. 
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This scaling behaviour was predicted in |12j using kinematic theory, and in pQ using 
equivariant bifurcation theory. The change from a rigidly rotating spiral to a travelling 
wave finger was described as a so called drift bifurcatior0 which occurs in the group 
dynamics rather than in the shape dynamics. To understand the bifurcation from 
rigidly rotating spirals to retracting fingers using the symmetry reduction method, 
one needs to look at the assumptions necessary for the orthogonal splitting of the 
full dynamics into the shape dynamics and the group dynamics - which underlies 
the freezing method as well as the theory in [T|. Symmetry reduction relies on the 
existence of a centre manifold. Finite-dimensional centre manifold reductions can be 
proven for spiral waves in unbounded domains assuming the existence of a spectral gap 
[37]. The spectral gap corresponds to a non-zero distance of the essential spectrum, 
which consists of the complement of the spectrum of the set of isolated eigenvalues of 
finite multiplicity, and the imaginary axis, thereby assuring normal hyperbolicity of 
the solutions. 

In [38] it was shown that this gap, in fact, does not exist for rigidly rotating 
spirals in unbounded domains. It was shown, however, that small perturbations to 
the unboundedncss of the domain, i.e. spiral wave solutions in 'very large' domains 
with imposed boundary conditions, open up a large spectral gap. This is in contrast 
to the situation for roll solutions say, where the spectral gap obtained in this way 
is negligible. Hence for many purposes it is reasonable to proceed as if a spectral 
gap is present so that a symmetry reduction can be performed. Such an approach 
has proved useful in understanding the transition to meandering and linearly drifting 
spirals [3j [51] [14j [TjJ [36] [3JJ [38] . This gap, however, becomes smaller the closer one is 
to criticality, at which point the spectral gap becomes zero. Close to criticality, when 
the overall curvature of the spiral wave becomes zero, and the solution has morphed 
into a semi-infinite travelling finger, finite-dimensional centre manifold theory is not 
applicable anymore and ought to be replaced by an infinite-dimensional Ginzburg- 
Landau type description. A rigidly rotating spiral wave may become unstable to an 
infinite number of modes of the continuous spectrum which cannot be captured by the 
freezing method or the bifurcation theory of pQ. However, we argue that the (possibly 
unstable) solution obtained by the finite dimensional reduction will, at least, function 
as an organizing centre for the full dynamics which takes into account the interactions 
with the continuous spectrum. 

Close to criticality, the break-down of the finite-dimensional description manifests 
itself in the requirement for an ever increasing resolution and accuracy in numerical 
simulations. Phenomenologically, one needs to resolve greater and greater parts of 
the spiral wave close to criticality, to resolve the behaviour of the far field spiral 
wave. A too small segment of the solution appears like a travelling wave without 
curvature. The tip region (i.e. the region which is not described by simple geometric 
constructs such as Archimedean spirals or involutes of circles) grows when criticality 
is approached, as was already encountered for small core spirals (cf. Figure 14.31 and 
Figure 15. 5|) . 

We found experimentally, that the actual values of r c and w at criticality, as 
well as the actual critical value of the excitability e c depend strongly on the numerical 
resolution. For example, the orientation of the finger within the computational domain 
has a strong influence, especially for moderate sizes of the computational domain. For 
a sufficiently small length of a spiral wave solution - which then appears finger-like, 



Note that the term drift bifurcation is also used in a different context for a pitchfork bifurcation 
of stationary patterns with reflection symmetry in an 0(2)-system 1281 1271 . 
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see Figure 14.51 - the proportion of the region dominated by the inaccurate Neumann 
boundary conditions becomes unproportionally large. We therefore rotate the finger 
in a pre-processing procedure to maximize the validity of the NBC. Spiral waves 
which do not leave the computational domain perpendicularly, are rotated around the 
centre of the rectangular computational domain, and then subsequently mapped back 
onto the computational grid using bilinear interpolation before the application of the 
freezing procedure. This procedure inevitably leaves triangles of the computational 
domain which have not been assigned values for the fields u and v after the rotation. 
We manually set u and v to zero on grid points falling into those triangles which 
is justified in the large core limit for small angels of rotation. Close to criticality 
it proved useful to use an iterative procedure, whereby reorienting and the freezing 
procedure are alternated at fixed e. However, this method of reorientation of spiral 
waves is inaccurate and not methodological, and small changes will have measurable 
effects in the values of the group parameters [53] . 

Similarly, the numerical results become more sensitive to the actual length of 
the finger which is resolved within the computational domain. At criticality, infinite 
resolution is required. This is illustrated in Figure 15.91 where we show results of the 
core radius function of e for simulations differing only in the length of the 

resolved finger within the domain. Whereas the values of r c are independent on the 
resolved size away from the bifurcation (but note the already large magnitude of r c 
in the log-scale), they differ strongly approaching criticality. The critical value e c also 
depends on the actual length of the resolved finger. This sensitivity of the results 
to the resolution is inherent and cannot be avoided due to the breakdown of the 
assumptions underlying the freezing method. 

It is pertinent to mention that, despite the sensitivity of the actual numerical val- 
ues of the rotation frequency, the core radius and the critical excitability, the linear 
scaling regime depicted in Figure 15.81 is robust against changes of size of the compu- 
tational domain, the orientation of the spiral and the discretization. 

The freezing method finds frozen solutions beyond the critical e c . In Figure 15.101 we 
show results for the core radius and the rotation frequency when the excitability is 
varied to values e > e c . We find that the linear scaling regime extends past the critical 
value e c . Moreover, the rotation frequency becomes negative for e > e c . Note that the 
theory of [T] does not describe the behaviour of solutions past the bifurcation point. 
For e > 0.0835 we were not able to find frozen solutions of the stationary system 
(|3.28j) . We have checked that the value of e = 0.0835 corresponds to the saddle node 
of travelling waves in ID, and exhibits the typical square root behaviour close to the 
saddle node bifurcation. 

At the bifurcation to retracting fingers at e c there is no guarantee that the so- 
lutions of the freezing method correspond to actual solutions of the original Barkley 
model. However, we have verified that the frozen solutions obtained for e > e c corre- 
spond to the retracting fingers by using them as initial conditions in the full Barkley 
model (|2Tjl . 

6. Summary. Our aim in this work was two-fold. In a first part, we have for- 
mulated a modification of the freezing method, introduced in . We have formulated 
the freezing method in polar and Cartesian coordinates for the time-dependent and 
the stationary formulation to freeze spiral waves in excitable media with non-diffusive 
inhibitors, typical for applications in cardiac dynamics. 

In particular, we have proposed a simple method to overcome oscillations near 

28 







24 



M6 



O 



-8 



o 



-16 



-25 -20 -15 -10 -5 

log(e - e c ) 



-24[ r .--- 

-25 -20 -15 -10 -5 

log(e - e c ) 



Fig. 5.8. Scaling behaviour of core radius r c (left) and rotation frequency u) (right). We show 
results for frozen solutions of the stationary system \3. 28V using Cartesian coordinates and NBCs, 



in computational domains with L x 
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e c = 0.08054091. 



The critical value of the excitability parameter is 
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Fig. 5.9. Demonstration of sensitivity of the critical excitability on the length of the resolved 
finger solution. Results are obtained by freezing the same spiral wave on rectangular domains with 
L x = 50 and L y = 62.5 (•), L y = 80 (x), L y = 120 (o) and L y = 140 (o). The resolved finger 
lengths are approximately Lq2.5 ~ 39.5 (•), Lgo f» 59 (x ), L120 ~ 97 (o) and L140 ~ 115 (o). 



the boundary, which have so far obstructed the investigation of the large core limit. 
Oscillations in the interior of the computational domain in the time-dependent prob- 
lem were eliminated by employing a semi-implicit Crank-Nicolson scheme. We have 
further introduced spiral boundary conditions by using Archimedean spirals and in- 
volutes of circles as geometric constructs to approximate contour lines of spiral wave 
solutions in unbounded domains and in computational domains whose size is much 
smaller then the actual physical domain. 

We have established the regions of applicability of our method. We found that to 
study spiral waves in the small core limit polar coordinates with SBCs are favourable, 
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Fig. 5.10. Core radius r c (left) and rotation frequency u> (right) for values of the excitability 
e > e c . At criticality the linear scaling behaviour is assumed from both sides of the bifurcation 
point e c = 0.08054091. We show results for frozen solutions of the stationary problem i3.28\ ) using 
Cartesian coordinates and NBCs, in computational domains with L x = 50 and L y = 62.5. 



whereas to study spiral waves in the large core limit, Cartesian coordinates and NBCs 
should be used. 

In a second part of this work, we have numerically investigated the large core 
limit of spiral waves. We have determined the shape of solutions near criticality, and 
have determined their rotation frequency as well as their core radius. Further, we 
discussed solutions of the freezing method for excitabilitics beyond criticality, which 
could be extended to the saddle node bifurcation of travelling waves. We have pre- 
sented results on the scaling behaviour of the spiral wave parameters in the large core 
limit and confirmed the linear scaling at the drift bifurcation developed in PQ, and 
discussed the limitations of the freezing method. We believe that these results may 
be helpful in designing kinematic theories. 

After submission of this work we became aware of work |16j in which the large core 
limit is investigated using a similar numerical method. The authors also identify a 
linear scaling regime in the large core limit, and study further meandering spirals 
and electrophoresis. In this work a phase condition is used which pins the tip of the 
spiral wave to the centre of the domain. We believe that this phase conditions with 
our implementation of the boundary conditions will prove useful in further studies of 
spiral wave dynamics. 
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